Allows you to define a complex (bioinformatic) workflow as a series of steps that involve input files and output files.
It identifies the dependencies between the steps and then runs all the steps needed to create a requested output file.
This greatly simplifies the orchestration of bioinformatics, and makes it much easier to find and re-run failed jobs.
Incredibly valuable for reproducible research:
Not just so others can reproduce your results
Also useful for you to quickly run your workflow on different clusters, etc.
That sounds pretty jargony!
Illustrate with an example
Hope that it piques the curiosity of some
Our Small Example: GATK Best Practices “Light”
flowchart TD
A(fastq files from 3 samples: our raw data) --> B(Trim the reads: fastp)
B --> C(Map the reads to a reference genome: bwa-mem2)
C --> D(Mark PCR and optical duplicates: gatk MarkDuplicates)
D --> E(Make gVCF files for each sample/chromo: gatk HaplotypeCaller)
E --> F(Load gVCFs into Genomic DB for each chromo: gatk GenomicsDBImport)
F --> G(Create VCFs from Genomic DB for each chromo: gatk GenotypeGVCFs)
G --> H(Concatenate chromosome-vcfs into a single vcf: bcftools)
A mini data set that only takes about 5 minutes to run through the major steps of a GATK-like variant calling workflow
Chinook salmon sequencing reads (a subset of our course example data).
Three paired-end fastqs from samples A, B, and C and data only from four chromosomes.
We will trim it, map it, mark duplicates, then make one gVCF file for each combination of individual and chromosome (only two chromosomes).
Then, call variants on each of two chromosomes.
Then catenate the resulting VCFs into a single VCF file.
Setting up our workspaces
Sync your fork’s main branch then pull updates into the main branch of your local clone.
Then, we simply need to ensure that we have a recent version of Snakemake.
We will create a mamba/conda environment snakemake-8.5.3 that has the latest version of snakemake as of this writing.
Once created, we activate that environment.
cd to the Snakemake-Example data directory inside the con-gen-csu repo.
On Alpine, get 4 cores on acompile.
# Create snakemake envmamba create -n snakemake-8.5.3 -c conda-forge -c bioconda snakemake==8.5.3# activate envconda activate snakemake-8.5.3# To make sure snakemake is working, print the help information# for snakemakesnakemake--help# If you want to make DAGs and rulegraphs, add GraphViz to this envmamba install anaconda::graphviz# change directories into Snakemake-Examplecd Snakemake-Example/# get onto acompile for interactive workacompile-n 4
Initial Configuration of our work directory
We can use the Unix tree utility to see what the Snakemake-Example directory contains.
Within the Snakemake-Example directory, type tree at the command line. This shows:
A Snakefile. Much more about that later.
A directory data with three pairs of FASTQ files
A directory envs that has information to install necessary software with conda
A directory resources that has genome.fasta: a FASTA file with the reference genome
flowchart TD
H(fastq files from 3 samples: our raw data) --> I(Trim the reads: fastp)
I --> J(Map the reads to a reference genome: bwa-mem2)
Some pseudo-shell code
# cycle over fastqs and do the trimmingfor S in A B C;dofastp-i data/${S}_R1.fastq.gz -I data/S{S}_R2.fastq.gz \-o trimmed/${S}_R1.fastq.gz -O trimmed/${S}_R2.fastq.gz \ other-arguments-etc...done# cycle over trimmed fastqs and do the mappingfor S in A B C;dobwa-mem2 mem resources/genome.fasta \ trimmed/${S}_R1.fastq.gz trimmed/${S}_R2.fastq.gzdone
What are some issues here?
Ah crap! I forgot to index genome.fasta!
This does not run the jobs in parallel!
Possible solutions for #2?
You can get things done in parallel using SLURM’s sbatch (which you probably need to use anyway).
Going about doing this with SLURM (a sketch…)
Consider the first two “steps”
flowchart TD
H(fastq files from 3 samples: our raw data) --> I(Trim the reads: trimmomatic)
I --> J(Map the reads to a reference genome: bwa mem)
Some pseudo-shell code
# cycle over fastqs and dispatch each trimming job to SLURMfor S in A B C;dosbatch my-trim-script.sh $Sdone# ONCE ALL THE TRIMMING IS DONE...# cycle over trimmed fastqs and dispatch each mapping job to SLURMfor S in A B C;dosbatch my-map-script $Sdone
What is not-so-great about this?
I have to wait for all the jobs of each step to finish
I have to explicitly start each “next” step.
If some jobs of a step fail, it is a PITA to go back and figure out which ones failed.
The dependence between the outputs of the trimming step and the mapping step are implicit based on file paths buried in the scripts, rather than explicit.
The Advantages of Snakemake
The dependence between input and output files is explicit
This lets snakemake identify every single job that must be run—and the order they must be run in—for the entire workflow (all the steps)
This maximizes the number of jobs that can be run at once.
The necessary steps are determined by starting from the ultimate outputs that are desired or requested…
…then working backward through the dependencies to identify which jobs must be run to eventually get the ultimate output.
This greatly simplifies the problem of re-running any jobs that might have failed for reasons “known only to the cluster.”
Snakemake is a program that interprets a set of rules stored in a Snakefile
Some explanations:
Rule blocks: the fundamental unit
Correspond to “steps” in the workflow
Keyword “rule” + name + colon
Indenting like Python/YAML
Typically includes sub-blocks of input, output, and shell
(Screen grab from Sublime Text which has great highlighting for Snakemake)
The rule:
Requires the input file resources/genome.fasta
Produces the output file resources/genome.dict
Uses software specified in envs/bwa2sam.yaml
Writes to a log file in results/logs/genome_dict.log
Executes the shell code samtools dict {input} > {output} 2> {log} to get the job done
What are those purple bits? {input}, {output}, and {log}?! in the shell code?
That is the syntax snakemake uses to substitute the values in the output, input, or log blocks (or other blocks…) into the Unix shell command.
Big Note: Output and log information is not written automatically to the output file and log file, nor is input taken automatically from the input file—you have to dicate that behavior by what you write in the shell block!
Thus, when this rule runs, the shell command executed will be:
We “drive” Snakemake by requesting the creation of output files
These output files are sometimes referred to as “targets”
snakemake looks for and uses the Snakefile in the current working directory.
Option -n tells snakemake to do a “dry-run:” (Just say what you would do, but don’t do it!)
Option -p tells snakemake to print the shell commands of the rules.
Those two options can be combined: -np
And we request resources/genome.dict as a target by just putting it on the command line:
Paste this into your shell
snakemake-np resources/genome.dict
And the output you got from that should look like:
What the output should look like
Building DAG of jobs...Job stats:job count------------------genome_dict 1total 1Execute 1 jobs...[Mon Feb 26 12:01:32 2024]localrule genome_dict:input: resources/genome.fastaoutput: resources/genome.dictlog: results/logs/genome_dict.logjobid: 0reason: Missing output files: resources/genome.dictresources: tmpdir=<TBD>samtools dict resources/genome.fasta > resources/genome.dict 2> results/logs/genome_dict.logJob stats:job count------------------genome_dict 1total 1Reasons:(check individual jobs above for details)missing output files:genome_dictThis was a dry-run (flag-n). The order of jobs does not reflect the order of execution.
Direct snakemake to create all the conda environments to run this workflow
We are going to start by telling snakemake to create the conda environments needed to run the whole workflow:
The results/vcf/all.vcf.gz is the ultimate target of the workflow and the --conda-create-envs-only tells snakemake to do nothing more than create all the conda environments needed to create the output file results/vcf/all.vcf.gz.
The output snakemake gives you should look something like this:
Expected output of the above
Building DAG of jobs...Your conda installation is not configured to use strict channel priorities. This is however crucial for having robust and correct environments (for details, see https://conda-forge.org/docs/user/tipsandtricks.html). Please consider to configure strict priorities by executing 'conda config --set channel_priority strict'.Creating conda environment envs/bwa2sam.yaml...Downloading and installing remote packages.Environment for /Users/eriq/Documents/git-repos/con-gen-csu/Snakemake-Example/envs/bwa2sam.yaml created (location: .snakemake/conda/03ade215a4c713db728206723f154153_)Creating conda environment envs/gatk.yaml...Downloading and installing remote packages.Environment for /Users/eriq/Documents/git-repos/con-gen-csu/Snakemake-Example/envs/gatk.yaml created (location: .snakemake/conda/d5b5e2cc0497b65f32514e037bc26ef9_)Creating conda environment envs/bcftools.yaml...Downloading and installing remote packages.Environment for /Users/eriq/Documents/git-repos/con-gen-csu/Snakemake-Example/envs/bcftools.yaml created (location: .snakemake/conda/cbf79ac1a639509f529e96f71cf3c38b_)Creating conda environment envs/fastp.yaml...Downloading and installing remote packages.Environment for /Users/eriq/Documents/git-repos/con-gen-csu/Snakemake-Example/envs/fastp.yaml created (location: .snakemake/conda/e37739b7128de56967cf47faf4c5cfef_)
Where do those environments get created? They do not go into your main conda library. Rather, as the output above shows, each environment gets stored inside .snakemake/conda in the current working directory.
A note on specifying conda environments for each rule
Before it runs the genome_dict rule, Snakemake will make sure that a conda environment with bwa-mem2 and samtools is installed in the .snakemake/conda directory. If it is not, then it installs it.
Then it runs the rule. Note that software installation happens only once.
Now, let’s run the genome_dict rule!
Go back and do the dry-run for the genome_dict rule again:
Paste this into your shell
snakemake-np resources/genome.dict
Now, to do a real run (not a dry run) we remove the -n (dry-run) option, and, because we are having snakemake manage the needed software using conda, we add the --use-conda option, and we tell it to use 1 core for those run: --cores 1. So, our command now looks like:
The output you get looks like what you saw before, but in this case the requested output file has been created.
And a log capturing stderr (if any) was created:
Paste this into your shell to see all the files
tree .
The output shows those two new files that were created
Output should look like this:
.├── Snakefile├── data│ ├── A_R1.fastq.gz│ ├── A_R2.fastq.gz│ ├── B_R1.fastq.gz│ ├── B_R2.fastq.gz│ ├── C_R1.fastq.gz│ └── C_R2.fastq.gz├── envs│ ├── bcftools.yaml│ ├── bwa2sam.yaml│ ├── fastp.yaml│ └── gatk.yaml├── resources│ ├── genome.dict <--- THIS IS A NEW FILE│ └── genome.fasta└── results└── logs└── genome_dict.log <--- THIS IS A NEW FILE6 directories, 14 files
Once a target file is created or updated Snakemake knows it
If you request the file resources/genome.dict from Snakemake now, it tells you that the file is there and does not need updating.
Paste this into your shell
snakemake-np resources/genome.dict
Because resources/genome.dict already exists (and none of its dependencies have been updated since it was created, and the code that creates the file in the rule has not been modified) Snakemake tells you this:
Expected output from Snakemake
Building DAG of jobs...Nothing to be done (all requested files are present and up to date).
This helps you to not remake output files that don’t need remaking!
Wildcards: How Snakemake manages replication
Wildcards allow running multiple instances of the same rule on different input files by simple pattern matching
If we request from Snakemake the file results/trimmed/A_R1.fastq.gz,
then, Snakemake recognizes that this matches the output of rule trim_reads with the wildcard {sample} replaced by A.
And Snakemake propagates the value A of the wildcard {sample} to the input block.
Thus Snakemake knows that to create results/trimmed/A_R1.fastq.gz
it needs the input files:
data/A_R1.fastq.gz
data/A_R2.fastq.gz
Try requesting those trimmed fastq files
See what snakemake would do when you ask for results/trimmed/A_R1.fastq.gz.
Paste this into your shell
snakemake-np results/trimmed/A_R1.fastq.gz
Note that you can request files from more than one sample:
Note that it will go ahead and start all those jobs independently, and concurrently, because they do not depend on one another. This is how Snakemake manages and maximizes parallelism.
Important Notes I: Multiple wildcards can be used together
Here, chromo and sample are two different wildcards. And snakemake understands that a requested file with a path like:
results/gvcf/NC_037122.1f5t9/A.g.vcf.gz
matches the output file:
results/gvcf/{chromo}/{sample}.g.vcf.gz
with
chromo = NC_037122.1f5t9
sample = A
Important Notes II: expand()
Snakemake provides some functions that are useful for creating lists of input or output files. Note that, at the top of the Snakefile we define python lists SAMPLES and CHROMOS.
# run it over 3 samplesSAMPLES = ['A', 'B', 'C']# variant calling over two "chromosomes"CHROMOS = [ "NC_037122.1f5t9", "NC_037123.1f10t14"]
Note that the {s} in curly braces get expanded according to s=SAMPLES, and the {{chromo}} gets turned into {chromo}, so that it is interpreted as a wildcard in that file list. (In general you escape braces by doubling them.)
Important Notes III: multiext()
Snakemake provides the multiext() function which is like expand() but for file extensions.
If Snakemake does not find a required input file for a rule that provides a requested output, it searches through the outputs of all the other rules in the Snakefile to find a rule that might provide the required input file as one of its outputs.
It then schedules all the necessary rules to run.
This means that an entire workflow with thousands of jobs can be triggered by requesting a single output file.
Short Group Activity
Trace the rules needed if we request the file results/vcf/all.vcf.gz.
One person from each group, write down the rule dependencies.
Do this by finding the rule that creates results/vcf/all.vcf.gz as output, then finding the rules that would create the input for that rule, and so on and so forth, all the way back to the original fastq files.
If any inputs change, Snakemake will re-run the rules that depend on the new input
Imagine that the sequencing center calls us to say that there has been a terrible mistake and they are sending you new (and correct) versions of data for sample C: C_R1.fastq.gz and C_R2.fastq.gz
Snakemake uses file modification dates to check if any inputs have been updated after target outputs have been created.
So we can simulate new fastq files for sample C by using the touch command to update the fastq file modification dates:
Paste this into your shell
touch data/C_R1.fastq.gz data/C_R2.fastq.gz
Now, when we run Snakemake again, it tells us we have to run more jobs, but only the ones that depend on data from sample C. Do a dry run to check that:
Paste this into your shell
snakemake-np results/vcf/all.vcf.gz
Check that it will not re-run the trimming, mapping, and gvcf-making steps for samples A and B, which are already done.
Snakemake makes it very easy to re-run failed jobs
Clusters and computers fail (sometimes for no apparent reason) occasionally
If this happens in a large, traditionally managed (Unix script) workflow, finding and re-running the failures can be hard.
Example: 7 birds out of 192 fail on HaplotypeCaller because those jobs got sent to nodes without AVX acceleration.
Five years ago, setting up custom scripts to re-run just those 7 birds could cost me an hour—about as much time as it takes me now to set up an entire workflow with Snakemake.
On the next slide we are going to create a job failure to see how easy it is to re-run jobs that failed with Snakemake.
Simulating a job failure as an example
First, let’s remove the entire results directory, so that we have to re-run most of our workflow.
Paste this into your shell
rm-rf results
Now, we are going to corrupt the read-2 fastq file for sample A (but keeping a copy of the original)
Now, run it with 4 cores and give it the --keep-going command which means that even if an error occurs on one job, all the other jobs that don’t depend on outputs from the failed job will still get run.
Snakemake runs as far as it can and then wraps it up, telling us that 8 of the 14 jobs were successful but at least one job failed:
Snakemake's concluding comments:
10 of 20 steps (50%)doneExiting because a job execution failed. Look above for error messageComplete log: .snakemake/log/2024-02-26T161639.704830.snakemake.logWorkflowError:At least one job did not complete successfully.
Here is a related, fun tip: Snakemake writes the log of every real (i.e., without the -n option) run into a log file in .snakemake/log. Try this:
ls .snakemake/log
If you want to get the log from the most recent run you can throw down some Unix:
You can use that to find any lines that say Error in rule in them (and 10 or 11 lines after it says “Error in rule”) to tell you more about the job that failed. e.g.:
grep"Error in rule"-A 11 $(ls-l .snakemake/log/*.log |tail-n 1 |awk'{print $NF}')
Re-running failed jobs is as simple as just re-starting Snakemake
First off, after a failure like that, we can always immediately do a dry-run to see what Snakemake must still do to finish out the workflow:
Paste this into your shell
snakemake-np results/vcf/all.vcf.gz
It is only going to require 10 more jobs to produce results/vcf/all.vcf:
Makes a graph (DAG) of the files involved in the workflow. If you view it, it looks like the following, which is exactly the sort of figure you would have created had you been extremely diligent about the tracing all the dependencies in the Snakefile with you team:
This tells you how many seconds it ran (s), the max RAM use (max_rss), amount of data read from disk and written to disk, total amount of CPU time, etc.
Super helpful
Easy to extract
Benchmark Example
Jobs processing 43 rockfish on NMFS on-premises cluster (SEDNA) vs in the cloud on Microsoft Azure (AZHOP)
Benchmark Example
Super easy to extract that for all the rules with a little Unix and R
This provides a complete BWA-GATK workflow including an arbitrary number of “bootstrapped-BQSR” rounds (which turn out to be pretty useless…)
Letting Snakemake interface with SLURM
The way to most easily allow Snakemake to dispatch jobs via the SLURM scheduler is by way of the Snakemake cluster option provided in a Snakemake profile.
A Snakemake profile is a YAML file in which you can record command line options (and their arguments) for Snakemake.
There is an officially supported Snakemake profile for SLURM, but I am partial to the Unix-based (as opposed to Python-based) approach to SLURM profiles for Snakemake described at: https://github.com/jdblischak/smk-simple-slurm.
A Snakemake profile for SLURM on Alpine
Alpine is the new NSF-funded cluster in Colorado
A Snakemake profile is simply a collection of command line arguments stored in a YAML file.
The set-resources and threads blocks are specific to my lcWGS workflow
Contents of hpcc-profiles/slurm/alpine/config.yaml